arXiv:1502.01562v2 [physics.flu-dyn] 4 Aug 2015 


Development of high vorticity structures in incompressible 3D 

Euler equations 

D.S. Agafontsev^’^, E.A. Kuznetsov^’^’^ and A.A. Mailybaev^’® 

^ P. P. Shirshov Institute of Oceanology, 

36 Nakhimovsky prosp., 117218 Moscow, Russia 
^ Novosibirsk State University, 2 Pirogova, 630090 Novosibirsk, Russia 
^ P.N. Lebedev Physical Institute, 53 Leninsky Ave., 119991 Moscow, Russia 
^ L.D. Landau Institute for Theoretical Physics, 

2 Kosygin str., 119334 Moscow, Russia 

^ Instituto Nacional de Matemdtica Pura e Aplicada - IMPA, Rio de Janeiro, Brazil 
^ Institute of Mechanics, Lomonosov Moscow State University, Russia 

We perform the systematic numerical study of high vorticity structures that de¬ 
velop in the 3D incompressible Euler equations from generic large-scale initial con¬ 
ditions. We observe that a multitude of high vorticity structures appear in the 
form of thin vorticity sheets (pancakes). Our analysis reveals the self-similarity of 
the pancakes evolution, which is governed by two different exponents and 

U/K describing compression in the transverse direction and the vorticity growth 
respectively, with the universal ratio ps 2/3. We relate development of these 

structures to the gradual formation of the Kolmogorov energy spectrum E]^ oc 
which we observe in a fully inviscid system. With the spectral analysis we demon¬ 
strate that the energy transfer to small scales is performed through the pancake 
structures, which accumulate in the Kolmogorov interval of scales and evolve ac¬ 
cording to the scaling law cjmax for the local vorticity maximums cjmax and 

the transverse pancake scales 1. 


I. INTRODUCTION 

The problem of whether incompressible 3D Euler equations develop a singularity in finite 
time (usually termed as blowup or collapse) from smooth initial data of finite energy is one of 
the most long-standing open questions in fluid dynamics and applied mathematics [TH5] . In 
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physics, formation of a collapse is considered as the most effective mechanism enhancing the 
energy dissipation of regular motion. For example, dissipation of oceanic surface waves takes 
significantly less time than its estimate due to viscosity [BHH]- The difference is attributed 
to formation of white caps on wave crests and snbseqnent breaking of waves, which is the 
manifestation of collapse for surface waves. For developed hydrodynamic turbulence, this 
process is connected with the energy transfer from large (forced) scales to small (viscous) 
scales, where the energy eventually dissipates. In the Kolmogorov-Obukhov theory puni, 
the velocity fluctuations 5v at intermediate spatial scales £ obey the power-law (|(^u|) oc 
where e is the mean energy flux from large to small scales. This formula can be 
easily obtained from the dimensional analysis, see e.g. CD [12]. According to the same 
dimensional argnments, the fluctuations of vorticity field uj = rot v diverge at small scales 
£ ^ 0 as (|(5a;|) oc while the time of the energy transfer from the energy-contained 

scale L to the viscous ones is finite and estimated as T ^ . Thus, formation of 

the Kolmogorov spectrum may be regarded to collapses of vorticity field; such mechanism 
can be observed in simplified (shell) models of turbulence [HI [Tl|. However, finite-time 
singularities are not necessary for the Kolmogorov tnrbulent spectra, since in case of the 
finite Reynolds numbers the exponential growth of vorticity is snfhcient for energy to reach 
viscons scales in finite time na. See also [iniini for nnmerical studies on early stages of 
turbulent spectra formation. 

There are several blowup and no-blowup criteria for the inviscid flows, which are useful 
both in analytical studies and numerical simulations. The widely used criterion is due to the 
Beale-Kato-Majda theorem [T8|, which states that the time integral of maximum vorticity 
must explode at a singular point. Several criteria, which also use the direction of vorticity, 
are developed by Constantin et al. ng, Deng et al. [201E] and Chae See also [221- 
I2H] for other regularity criteria. The blowup scenario based on the vortex lines breaking 
(or overturning) was analyzed in [20] in the framework of the integrable incompressible 
hydrodynamic model with the Hamiltonian J \cu\dr and the same symplectic operator [22] 
as for the 3D Enler equations (such unusual Hamiltonian can be obtained from the 3D Euler 
equations in the so-called local indnction approximation). Eormation of a singularity in the 
framework of renormalization group formalism was discussed in [2EH2D • 

Possible formation of a singularity in incompressible 3D Euler equations has been exten¬ 
sively studied in the past decades with direct numerical simulations, which mainly refer to 
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the flow in a box with periodic boundary conditions. Below we provide a short review of 
these numerical results; see also a brief but extensive account in [32], as well as [331 El] for 
the studies of blowup triggered by the boundary. 

In 1992, Brachet et al. [35] using spectral methods studied periodic flows with random 
initial conditions on 256^ grids and also the Taylor-Green vortex on 864^ grids, and found 
the energy spectrum well-approximated as Ek{t) = . The exponent 8{t) clearly 

demonstrated the exponential decay with time, 5{t) oc Maximum of vorticity was 

growing almost exponentially in time with the total increase by the factor of 6 for both 
initial conditions. The regions of high vorticity represented pancake-like structures in phys¬ 
ical space, which were compressed in one direction while their sizes in other two directions 
did not change considerably. Such tendency towards the vortex sheets leads to depletion 
of nonlinearity and prevents the formation of a finite-time singularity, see the related dis¬ 
cussion in [331133] . Thus, further numerical studies were mainly focused on specific initial 
conditions providing enhanced vorticity growth [33], e.g., antiparallel or orthogonal vortices. 
Initial conditions were usually chosen to be symmetric, as this required less computational 
resources. 

In 1993, Kerr [39] analyzed the interaction of two perturbed antiparallel vortex tubes on 
grids of up to 512 x 256 x 128 using the Chebyshev method. The results were interpreted 
in favor of the blowup, max|a;| ~ (to — with an increase of the vorticity maximum 
by the factor of 24. In the subsequent publication [10], the two characteristic length scales 
of the singularity were identified as p ~ (to — t) and R ^ {to — t)^/^. However, recent 
numerical simulations performed by Hou and Li [ 2 ll HI] with the pseudo-spectral method 
on the 1536 x 1024 x 3072 grid questioned the blowup behavior for these initial conditions. 
The solution was extended beyond the earlier estimated blowup time, showing that the 
maximum vorticity evolution is slower than doubly exponential. Later the analysis was 
reconsidered by Bustamante and Kerr [32] , suggesting the hypothesis of the vorticity growth 
as max I u; I ~ (to — t)”'*' for 7 > 1, and then by Kerr [33] concluding with the double¬ 
exponential growth. 

In 1998, Grauer et al. mi used the adaptive mesh refinement technique claimed to be 
equivalent to non-adaptive 2048^ grids. The authors achieved maximum vorticity increase 
of about 10 times and interpreted their results in favor of the blowup hypothesis, max |a;| ~ 
(to — t)~^. Similar conclusions were drawn in 2012 by Orlandi et al. [3H] for two colliding 
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Lamb dipoles on grids up to 1536^. However, in both cases m and [IS] the initial conditions 
had a singularity of vorticity derivatives. 

Several studies considered the Pelz-Kida initial flow of high symmetry with the indication 
of blowup behavior reported in earlier works In 2008, Grafke et al. [l9] performed 

comparison of different spectral and real-space numerical methods with the Pelz-Kida like 
initial flow. The authors achieved the increase of maximum vorticity by about 2.5 times 
for 512^ grids and about 3.5 times for 1024^ grids before the methods started to diverge 
noticeably. Then the adaptive mesh refinement simulations were carried out with the effec¬ 
tive resolution of 4096^ and vorticity maximum increase of about 6 times. These results, 
also confirmed by Hou and Li [50], demonstrated no tendency toward blowup at the times 
predicted earlier. It was also noted that the vorticity increased exponentially with time on 
the Lagrangian trajectory, which ended at the maximum vorticity point at the end of the 
simulations. 

In the present paper, we analyze numerically the 3D Euler equations for rather generic 
large scale initial conditions in a periodic box, i.e., focusing on typical development of 
vorticity in inviscid incompressible flows. For numerical simulations, we use the pseudo- 
spectral method on adaptive rectangular grid. Our goal is the systematic study of high 
vorticity structures, including the maximum vorticity region and other local phenomena. 
We identify these structures by looking at local maximums of the vorticity modulus. We 
show that a multitude of high vorticity structures appear and have the form of pancakes, i.e., 
thin vorticity sheets compressing exponentially in time in a self-similar way for the pancake 
transverse direction. Unlike the pancake model proposed in |35|, this self-similar dependence 
is governed by two different exponents for the pancake compression and the vorticity growth; 
we suggest how the pancake model [2H| can be modified to capture the observed behavior. 
Due to the exponential vorticity growth, no tendency toward the blowup is observed. 

Our main result is the demonstration of close relation between the collective evolution 
of pancake vorticity structures and the formation of Kolmogorov turbulent spectra. Our 
simulations show that the Kolmogorov energy spectrum starts to form in the fully inviscid 
system, where numerical error is kept very small (the flow is not affected by small-scale 
finite grid effects as opposed to Ell). With the analysis of local maximums we argue that 
the Kolmogorov spectrum is attributed to the pancake structures, which accumulate in 
the same interval of scales and evolve according to the scaling law for the local vorticity 
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maximums as cj^ax oc where £ is the pancake thickness. Though our conclusions are 

limited by the numerical resources that allow the Kolmogorov interval of about one decade, 
our results provide a new insight on the development of high vorticity structures in relation 
with the turbulent spectra, and indicate the importance of further study in this direction. In 
addition to these results, in our next paper we will present numerical simulations performed 
in the vortex lines representation variables. The vortex lines representation developed by 
Kuznetsov and Ruban [26l [52l |53] describes the flow from the point of view of the moving 
vortex lines, which are compressible. This representation helps in understanding some of the 
numerical observations presented in the current paper, including the scaling law (Umax oc 
between the local vorticity maximums (Umax and the transverse pancake scales £. 

The paper is organized as follows. Section |TI] describes the numerical method. In Sec¬ 


tion [in] we study the pancake structure near the global vorticity maximum, while in Sec¬ 
tion [I 3 we focus on the statistics of local maximums. Section m discusses the numerical 
observation of the Kolmogorov spectrum and its relation to extreme vorticity structures. 
The final Section |VI| contains conclusions. The Appendices [^ and [^ contain the initial 
conditions and the results of the second simulation. 


II. NUMERICAL METHOD 


The Euler equations describing dynamics of ideal incompressible fluid of unit density in 
three-dimensional space are 

dv 


dt 


+ (v • V)v = — Vp, divv = 0, 


( 1 ) 


where the velocity held v = {vx,Vy,Vz) and pressure p are smooth functions of space coordi¬ 
nates r = {x, y, z) and time t. We consider solutions of Eq. Q in the box r G [— vr, tt]^ with 
periodic boundary conditions. Eor numerical simulations, we use the formulation of Euler 
equations in terms of vorticity (also known as Helmholtz’s vorticity equations). 


dijj , . 

-^ = rot (v X cj), V = rot UJ. 
dt ^ 


( 2 ) 


Assuming the vanishing average velocity J vd^r = 0, the inverse of the rotor operator in 
Eq. 0 is uniquely defined and has the form 

ik X Cc;(k) 


v(k) = 




for k = ||k|| 7 ^ 0 ; v(0) = 0, 


(3) 
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FIG. 1: (Color on-line) Normalized spectra of the vorticity field Sj{k)/ max^?^ Sj{k) along the j-axis 
{x solid black, y dashed blue, z dash-dot red) for the initial data Ii (see Eq. Q and Tab. [^in the 
Appendix 0 8it t — 5 and at the final time t — 6.89 (inset). 

for the Fourier transformed velocity v(k) and vorticity cj(k). The wavevector k = {kx^ ky, k^) 
has integer components for the periodic box, r G [— 7r,7r]^. Note that the transformation (|^ 
automatically satisfies the incompressibility condition, k • v(k) = 0. 


A. Adaptive scheme and initial conditions 


We solve the system ([^ numerically using the pseudo-spectral method in the adaptive 
rectangular grid, which is uniform along each coordinate. The number of nodes N^, Ny and 
Nz in each direction is adapted independently, as explained below. To avoid the so-called 
bottle-neck instability we perform the filtering in Fourier space with the cut-off function 


p(k) = exp ( — 36 


kx 


Kx 


36 




k'li 


Kx 


36 




36n 


(4) 


Here Kj = Nj/2 are the maximum wavenumbers along directions j = x,y,z, i.e., kj G 
[—Kj,Kj]. Function (|^ cuts off approximately 20% of the spectrum at the edges of the 
spectral band in each direction. We use the fourth-order Runge-Kutta method with the 
adaptive time stepping, which is implemented according to the CFL stability criterion with 
the Courant number 0.5. 

For optimal performance, we start simulations with the 128^ grid, which is appropriate 
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for large-scale initial conditions. The adaptive scheme is controlled by the functions 

Sj{k) = j \u>{p)\^S{\pj\ - k)d^p, j = x,y,z, (5) 

which describe the enstrophy spectrum along each axis. As shown in Fig. these functions 
have breakpoints approximately at the level of ~ 10“^°, which corresponds to numerical 
noise. During the simulation, these breakpoints move to larger k reflecting the excitation of 
higher harmonics. We stop the simulation every time when one of the breakpoints reaches 
the value of 2Kj/3 for the corresponding axis, and then continue with the refined grid which 
has increased number of points Nj along the axis j. The vorticity field is interpolated to 
the new grid using the Fourier interpolation procedure, which has an error comparable to 
round-off. The simulations carried out in this way are not affected by the aliasing errors. 
After the total number of nodes N^NyNz reaches 1024^, we continue the simulation with 
the flxed grid and stop when any of the values Sj{2Kj/3) reaches 10“^^ max^ 5j(A:), see 
Fig. inset). 

We tested our simulations with smaller time steps, different harmonics filtering, includ¬ 
ing the standard 2/3 cut-off rule [H], and found negligible difference in the results. The 
simulations conserve the total energy E = | J |vp dr and helicity D = f (v ■ u;) dr with a 
relative error smaller than 10“^^. We also compared our simulations performed with differ¬ 
ent limitations for the maximal total number of nodes, and found that the results perfectly 
converge. For instance, in the time interval 4.13 < t < 5.83, where the simulations with 
N^NyNz < 512^ and N^NyNz < 1024^ run on different grids, the relative difference for the 
maximum vorticity is less than 10“^. The main source of this difference is attributed to 
the fact that different nodes of the grids are used for approximating the maximum. Note 
that the high accuracy used in our simulations is crucial for the detailed analysis of extreme 
structures in the flow, which requires very accurate computation of higher-order derivatives 
near local maximums of the vorticity. For less accurate simulations, numerical errors in 
these regions increase and the simulation results become unsatisfactory for our purposes. 

We consider initial conditions represented in the form of Fourier series 

t = 0 : tj(r) = [Ah cos(h • r) -|- Bh sin(h ■ r)] , (6) 

h 

where h = {hx,hy,hz) is a vector with integer components. Since divcj = 0, the real 
vectors Ah and Bh must satisfy the orthogonality conditions, h • Ah = h • Bh = 0. We 


fix the vectors A(o,o,i) = (0,1,0) and B(o,o,i) = (1,0,0), and choose the other coefficients as 
random numbers with zero mean and variance ^ exp(— |hp). Thus, the initial conditions 
represent the large-scale vorticity field given by the shear flow (degenerate ABC flow) 

LJx = sin(2), u)y = cos(2;), = 0, (7) 

which is the exact stationary solution of the 3D Euler equations, with a random perturbation. 

A number of initial conditions were tested with the purpose of choosing good candidates 
for the final high precision simulations. We selected the two initial conditions with better 
performance for the global vorticity maximum and relatively small number of excited har¬ 
monics, which are denoted as A and A and summarized in Tabs. |I] and [IT] in the Appendix [A| 
along with some simulation information. We focus our analysis on the simulation for the 
initial condition A, which reaches the time t = 6.89 on the grid 486 x 1024 x 2048. The 
different number of grid points for different directions results from the adaptive scheme, 
which resolves the anisotropy of vorticity in an optimal way. In the Sections below this 
simulation will be always assumed unless otherwise stated. The simulation results for the 
initial condition A are summarized in the Appendix 

B. Analysis of high vorticity structures 

As we demonstrate in the next Sections, the regions of high vorticity are strongly 
anisotropic and develop in the form of pancake-like structures, which are thin in one direction 
and remain large in other two directions. Fig. |^a). Conclusions of our paper rely signifi¬ 
cantly on the systematic identification of such structures, which is performed by searching 
for local maximums of the vorticity modulus |u?(r)|. The values of vorticity in the pancake 
are very sensitive to the distance from the pancake midplane, while the shifts along this 
plane have much smaller effects. Thus, searching for maximums by a simple comparison 
of neighboring grid points yields a large number (up to thousands) of “false” local maxi¬ 
mums. For this reason, we developed a three-step numerical procedure for the identification 
of “true” local maximums. 

At each grid point r^, we compute the gradient vector g(rg) = V|a>| and the Hessian 
matrix H(r 3 ) = [d‘^\u}\/dxidxj] using the spectral method. A local maximum is character¬ 
ized by the vanishing gradient vector and negative definite Hessian matrix, which has three 
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FIG. 2: (Color on-line) (a) Vorticity field |a;| at fixed x = xq at the final time t — 6.89. Here xq 
is the x-coordinate of the global maximum, (b) Positions (circles) of local maximums of |a;| and 
directions (lines) of the eigenvectors wi normal to the pancake structures, at t = 6.89. Red color 
marks the global maximum. 


negative eigenvalues Ai < A 2 < A 3 < 0 with the orthonormal eigenvectors Wi, W 2 , W 3 . As 
we will demonstrate below for the local maximums of vorticity, |Ai| IA 2 I IA 3 I and the 
eigenvector Wi determines the perpendicular direction to the pancake midplane while the 
eigenvectors W 2 and W 3 lie in this plane. Assuming that the local maximum lies near r^, 
the second-order approximation yields 

~r3-g(r3)H-i(r3). (8) 

As the first step in our procedure, we check every grid node to satisfy the following two 
conditions simultaneously: Eq. ([^ yields the point lying in one of the adjacent grid cells, 
and the Hessian matrix H(rg) is negative definite. These conditions determine a cloud of grid 
points around each local maximum. So, as the second step, we select a single node in each 
cloud by choosing the node with the smallest distance — r^l. At larger times some of the 
pancake structures contain several local maximums of vorticity. Some of these maximums 
form localized clusters, within which local maximums have close values of vorticity, almost 
the same eigenvectors Wi, and he approximately in the same plane almost perpendicular 
to the eigenvectors wi. Thus, such maximums represent close and similar oriented parts of 
the same pancake. In the third step we determine all such clusters of local maximums, and 
leave only one point with the largest value of vorticity from each cluster. 

We expect that the local maximums in the resulting set are associated with different 
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FIG. 3: (Color on-line) (a) Evolution of the global vorticity maximum (logarithmic vertical scale). 
The dashed line indicates the exponential slope oc with characteristic time — 2. (b) 

Evolution of characteristic spatial scales for the region of the global vorticity maximum, ii (black), 
(.2 (blue) and £3 (red). The dashed line indicates the exponential slope oc with characteristic 

time Ti = 1.4. 

pancake structures or different (distant or differently oriented) parts of the same pancake. 
Fig. [^b). The proposed method allows adequate selection of local maximums of vorticity, 
producing a little amount of false maximums. We compared its performance on significantly 
different grids, and found that the resulting sets of local maximums differ by no more than 
10-15%. 


III. EVOLUTION NEAR THE GLOBAL MAXIMUM OF VORTICITY 

Numerical studies of singularities in the 3 D incompressible Euler equations use the evo¬ 
lution of global vorticity maximum = max|u;(r)| as one of the principal tests. 

According to the Beale-Kato-Majda theorem [TH], unbounded increase of this maximum is 
necessary for the finite-time blowup. In this Section we consider the dynamics of ooraaxit) and 
analyze the local geometry of the flow. The global vorticity maximum grows with time from 
the initial value a;inax(0) ~ 1.5 to 11.8 at the finial time of the simulation t = 6.89, Fig. |^a). 
For t > 4.5, this growth is well-approximated by the exponential function oc 

with ~ 2 . Contrary to the vorticity, the velocity maximum max |v(r)| does not change 
more than by 10 % during the whole simulation. 

Local geometry near the vorticity maximum can be studied using the Hessian matrix 
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FIG. 4: (Color on-line) Isosurface of constant vorticity |u;| = O.Swmax in the local coordinates 
(ai, a 2 ) 0,3) at the final time of the simulation, t — 6.89. Note much smaller vertical scale. 

H = [d‘^\u:\/dxidxj] computed at the maximum point r^. This matrix has three negative 
eigenvalues Ai < A 2 < A 3 < 0 with the orthogonal eigenvectors wi, W 2 , W 3 . Considering 
the local orthonormal basis as r = + UiWi + a 2 W 2 + UsWs, the vorticity modulus can be 

described by the quadratic approximation 

^ / /I r. \ / /1o \ - I V.f, I _ 

(9) 


|a;(r)| 






£j = 


-A, 


The quantities £1, £2 and £3 determine the size of high vorticity region at fixed time, and 
their time evolution is shown in Fig. |^b). Note that the discontinuities in this figure near 
t = 2 correspond to the change of global maximum among different local maximums, while 
the noise in determining £2 and £3 for larger times is the result of amplified numerical error 
due to the ill-conditioned Hessian matrix. The smallest size decreases exponentially in time, 
£1 oc with 1.4 for t > 4.5, while £2 and £3 remain almost the same. At final 

time, we have £i /£2 £i /£3 ~ 10“^, which implies that the high-vorticity region represents 

a very thin pancake structure with the normal vector Wi and tangent vectors W2 and W3. 
This is confirmed in Fig. (note the much smaller scale of the vertical axis) showing the 
numerically computed isosurface |cj(r)| = O.StUmax in the coordinates (01,02,03). 

Similar pancake structures were observed systematically in j35], where the singularity 
was analyzed by looking at the analyticity strip of the solution. This method relies on the 
evolution of the energy spectrum 




( 10 ) 
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(b) 


(c) 


FIG. 5: (Color on-line) (a) Energy spectrum Ek{t) at different times. Figures (b) and (c) show the 
exponents n{t) and 5{t) of the fit Ek{t) « . The dashed fine in figure (c) indicates 

the exponential slope with characteristic time Ts — 1.4. 


where o is the spherical angle. Numerically, we find Ek{t) as a sum | ^ |v(k) p over all nodes 
in the spherical shell A: < |k| < fc + 1. This procedure is simple and yields the result, which 


is very close to the direct computation of the integral in Eq. (10) with the interpolation of 
velocity Fourier components on the sphere of radius k. 

For fixed time and sufficiently large wavenumbers, the energy spectrum decays with 
k exponentially, Ek{t) oc possibly with an algebraic prefactor, until the level of 

~ 10“^^ corresponding to numerical noise. Fig. [^a). The exponent S(t) can be associated 
with the width of analyticity strip for the solution extended to complex coordinates r. In 
[SS], the spectrum fitting Ek{t) ~ c{t)k~'^^*'> was used and the function 5{t) showed the 
exponential decay 8(t) oc at large times. Our simulations lead to the same conclusions 

with E ~ 1-4, Fig. 1^ c). Note that E ^ E, where E Ri 1.4 describes the exponential decay 
of the pancake width (we estimate the absolute numerical accuracy for the time scales 
E and E as ±0.1). This relation is natural, since the analyticity strip must be determined 
by the most extreme event, i.e., by the thinnest part of a pancake at the point of maximum 
vorticity. Note that our results for the algebraic prefactor n(t) differ from [33], where n(t) 
approached —4 at late times. Fig. [^b). 

The following self-similar solution of the Euler equations was suggested in [35] for the 
description of the flow structure in a small neighborhood of the pancake, 

2 _L ^2 

( 11 ) 


a\ 


a 2 „ / ,/rx a{ ± 

ni =^2 = —, V3 = f{aie^), p= - 
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which is written in local coordinates (ai, a 2 , as) with the axis ai perpendicular to the pancake; 
Vj are components of the velocity, p is pressure and /(ai) is an arbitrary function. The model 
0 leads to the exponential growth oc of vorticity uj = (wi, 6 <; 2 , <^ 3 ) as 

u}\ = = 0 , u)2 = ( 12 ) 


in the pancake, whose thickness decreases exponentially as oc with the same charac¬ 

teristic time T. However, our simulations demonstrate that the evolution of the pancake is 
governed by the two different exponents ~ 2, and T( ~ 1.4, for the vortic¬ 

ity growth and the pancake compression, respectively. We propose that the flow near the 
pancake can be approximated as 


vi = V 2 = ^, Us = e “V (oie*/^«) , p = 

If • - 


a? -h Un 

2Ti ^ 


(13) 


where a = l/T^ — l/T^ > 0. The modified model (13) is not the solution of the Euler 
equations. However, it satisfies the incompressibility condition divv = 0, yields the correct 
exponents for the maximum vorticity and the pancake thickness 


Ui=u}3 = 0, a;2 = -e*/^“/'(oie*/^^). 


(14) 


and after being substituted into the Euler equations ([^, ensures the cancellation of the 
leading-order terms growing exponentially as oc while the next-order terms (not can¬ 

celed) decay as oc e““h 

Our pancake self-similarity hypothesis agrees well with the simulation results, see Fig. 
Note that the vorticity vector near the pancake is mostly aligned with the eigenvector W 2 


(axis 02 ), in agreement with the model given by Eqs. (13)-(14). One can expect that the 


pancake (13)-(14) contributes to the energy spectrum in the interval of wavenumbers \k\ ^ ki 
with 

= W (15) 

where £1 oc is the characteristic scale along the axis ai in Eq. § and L estimates the 

oscillation period of vorticity components in Fig. |^a); this relation can be deduced, e.g., by 
analogy with the Fourier transform of the function sech(x/£i). At final time, the numerical 


simulation yields L 0.13 and ii ~ 0.018, which is in good agreement with relation (15). 


It is worth noting that self-similar model (13)-(14) cannot be interpreted in the context of 


two-dimensional flow, since vorticity grows exponentially in time, while the vorticity vector 
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FIG. 6 : (Color on-line) (a) Components of the vorticity vector m = wiwi + 0 ^ 2 W 2 + 0 ^ 3 W 3 along the 
axis ai perpendicular to the pancake &tt — 6.89. (b) Renormalized vorticity component LO^./co -m^^ cU) 
vs. renormalized coordinate ai/£i at different times demonstrating self-similarity of the vorticity 
field. 

is parallel to the pancake plane and is not orthogonal to the velocity. This indicates the 
importance of third dimension for the pancake dynamics. 

The second simnlation with the initial condition I 2 follows the same scenario for the 
global vorticity maximum and the associated region of high vorticity; see Appendix 
The corresponding characteristic times for the exponential behavior are estimated as ~ 
2.7, Ti ~ 1.7 and Ts ~ 1.7. Again, the two time scales Ts ~ Ti are close, while the 
scale is considerably larger. We conclude that the pancake behavior is governed by the 
two characteristic times controlling the exponential growth of vorticity and the exponential 
decrease of pancake thickness. 


IV. LOCAL MAXIMUMS OF VORTICITY 


Figure [^a) shows the vorticity distribution for the cross-section passing through the 
global maximum at final time. One can see that the regions of increased vorticity have a 
tendency to form a number of (thin and wide) pancake structures, which is in agreement with 
earlier simulations for generic initial conditions [HH] . We expect that the structural analysis of 
these pancakes might give insight into the 3D Euier dynamics. However, direct identification 
of the pancakes is a complicated numerical problem. We approach this problem by finding 
and analyzing local maximums of the vorticity modulus, as described in Section [II B[ 
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FIG. 7: (Color on-line) (a) Total number of local vorticity maximums versus time obtained for 
the two simulations corresponding to the total number of nodes NxNyNz < 512^ (black, ends at 
t — 5.83) and 1024^ (red), (b) Evolution of local vorticity maximums (logarithmic vertical scale). 
The dashed red line indicates the exponential slope oc with characteristic time = 2. 


Figure [^a) shows that the number of local maximums increases with time, from 3 at 
t = 0 to 16 at the end of the simulation. The values of vorticity modulus at maximum 
points are shown in Fig. [^b). The results suggest that the vorticity growth tends to be 
exponential, mmax(^) for many of the local maximums with rather close values of 

the characteristic times 

Figure [^a) shows the evolution of the three length scales, £i, £2 and £ 3 , computed for 
each local maximum according to Eq. 0 . We see that the smallest scale £\ (black) decays 
nearly exponentially, £1 oc with rather close values of the characteristic time Ti for 

different local maximums, while the scales £2 (blue) and £3 (red) remain the same or decrease 
slightly; recall that the large numerical error for the computation of £3 is related to the ill- 
conditioned Hessian matrix. This hgure conhrms the visual observation (Fig. that most 
of the regions of increased vorticity tend to have the pancake shape characterized by the 
small spatial scale (thickness) £1 exponentially decreasing with time. 

It is very instructive to see the distribution of local maximums among the spatial scales 
at hxed time, which highlights the contribution of different pancake structures to the en¬ 
ergy at different scales. Using the estimate (15), we plot in Fig. [^b) the characteristic 
wavenumbers ki for local maximums in decreasing order. The results suggest that the two 
leading maximums propagate much faster towards large wavenumbers. However, all the 
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FIG. 8 : (Color on-line) (a) Evolution of characteristic spatial scales (black), ^2 (blue) and ^3 
(red) of local maximums, see Eq. Exponentially decreasing values of indicate formation 
of the pancake structures in the vorticity field, (b) Characteristic wavenumbers ki estimated by 


Eq. (15) for the local maximums versus local maximum index numbers at different times. Local 
maximums are sorted in decreasing order of k\. Large circles mark the global maximums. 


other local maximums (i.e., most of the pancake structures) fill densely the interval from 
large to medinm scales. Thronghont this interval the distribution of maximnms has a well- 
established slope (i.e., the pancakes are distributed with a specific “spectral density”). This 
interval increases with time and reaches 0 < fci 30 at the final time t = 6.89. 


V. VORTICITY STRUCTURES AND THE KOLMOGOROV SPECTRUM 

In fully developed turbulent flow with large Reynolds numbers, the energy from large 
scales is transported through a wide range of medium scales (the inertial interval) to small 
scales, where it is eventually dissipated, see [fillT^. The viscosity is negligible in the inertial 
interval, where the evolution can be described by the Euler equations. The dimensional 
considerations suggest the well-known Kolmogorov scaling law, oc for the energy 

spectrum in the inertial interval. Our simulations can be seen as describing the initial stage 
of turbulent dynamics, developing from large-scale initial data, at times before the flow gets 
excited at viscous scales. Generally speaking, considerations of the Kolmogorov theory do 
not extend to this case, because the energy cascade is not formed yet. Thus, the scaling 
law Ek oc is not necessarily satisfied. From this point of view, our simulations of the 

3D Euler equations give the possibility for studying the initial stage in formation of the 
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FIG. 9: (Color on-line) Energy spectrum Ek{t) at different times. Straight line above the curves 
indicates the slope of the Kolmogorov power-law, Ek oc 

Kolmogorov spectrum together with the respective energy transfer mechanisms. 

Figure shows that, at sufficiently small wavenumbers, we clearly observe the gradual 
formation of the Kolmogorov interval Ek oc This interval grows with time and extends 

to a decade of wavenumbers, 2 20, at the end of the simulation. The Kolmogorov 

interval corresponds to the “frozen” part of the energy spectrum: Ek{t) changes slightly 
with time in the Kolmogorov region in contrast to the vast changes at larger wavenumbers. 
Taking into account the times and logarithmic scale in Figure]^ one can guess that the size 
of the Kolmogorov interval increases exponentially in time. 

In order to link the energy spectrum Ek with the pancake vorticity structures studied in 
the previous Sections, we look at the flow in Fourier space. One can expect that each pancake 
generates a structure extended in one direction in Fourier space (“jet”), aligned with the 
eigenvector wi perpendicular to the pancake (such jets form, e.g., in the two-dimensional 
inviscid flow |5ll |55]). Inside such a jet, the Fourier components of the flow should be large 
in comparison with the remaining background. In order to visualize this effect, we consider 
the function 

a;(k) = (^’(k) / max |u;(p)|, (16) 

/ |p|=fc 

representing the Fourier transformed vorticity oj (k) scaled to the maximal norm within each 
spherical shell. The reason for such a normalization is to compensate the strong decay of 


vorticity with k. Numerically, the maximum in Eq. (16) is computed among the nodes in a 
spherical shell of unit thickness. 
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FIG. 10: (Color on-line) Isosurface |u;(k)| = 0.2 of the normalized vorticity field ( [T^ in Fourier 
space at final time t — 6.89. (a) The jet is aligned with the eigenvector wi (dashed black line), which 
is the normal direction of the pancake structure at the global vorticity maximum in physical space, 
(b) The closer view with the eigenvector wi (solid black arrow) for the global vorticity maximum 
and the respective eigenvector (dash-dot blue arrow) for the third largest local maximum. 

Figure [To|(a) shows the isosurface |ci;(k)| = 0.2, where the very thin interior part cor¬ 
responds to larger vorticity, |d;(k)| > 0.2. As expected, this isosurface is aligned with the 
eigenvector wi computed at the global maximum, which should bring the dominant con¬ 
tribution to Fourier components of vorticity at large k. Figure lEb) demonstrates that 
the isosurface geometry is different at smaller wavenumbers, A: ^ 20, corresponding to the 
Kolmogorov interval in Fig. Here, different jets contribute and some of these jets can be 
clearly related to the pancakes of other local vorticity maximums (the figure shows directions 
for the first and third largest local maximums, while the second local maximum yields the 
direction very close to that for the first one). As shown in Fig. [TT|(a), the small interior part 
of this isosurface dominates in the energy spectrum. These observations reflect the extreme 
anisotropy of vorticity field and suggest that the Kolmogorov energy spectrum for /c < 20 
may be related to a collection of the pancake structures, rather than being determined solely 
by the dominant one. To further test this supposition, we integrate the energy spectrum in¬ 
side each of the two jets |th(k) | > 0.2 shown in Fig. [I0|(b), which can be separated from each 
other by cones of angle 0/2, drawn from the origin around the corresponding eigenvectors 
wi. Here 0 ~ 7r/4 is the angle between the directions of the jets. One can see from Fig.[TT|(b) 
that the contributions of the jets to the energy spectrum are comparable at sufficiently small 
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FIG. 11: (Color on-line) Energy spectrum at the end of the simulation, (a) Calculated inside 
the isosurface |ci;(k)| > 0.2 (solid black), outside this isosurface (dash-dot blue) and the sum of the 
two (dashed red), (b) Calculated inside the isosurface |ci;(k)| > 0.2 (thick solid black), inside the 
first jet (dashed blue), inside the second jet (thin solid green), and the sum of the two jets (dash-dot 
red). Black line above the curves indicates the slope of the Kolmogorov power-law, oc 

wavenumbers in the Kolmogorov region, while at larger wavenumbers (both inside and out¬ 
side the Kolmogorov region) the contribution from the leading jet becomes dominant. Thus, 
the energy spectrum at sufficiently large wavenumbers /c ^ 20 is determined mainly by the 
leading pancake structure corresponding to the global vorticity maximum. This means that 
the process of energy transfer to small scales is performed through the evolution of the pan¬ 
cake vorticity structures in the physical space. This interpretation is similar, in a spirit, 
to the simplified turbulence models [56] based, e.g., on Lundgren vortices [SZl EH]- Note, 
however, that our simulation describes the initial stage of turbulence formation, when basic 
assumptions of these models related to stationary developed turbulence do not apply. 

Statistical properties of the pancakes analyzed in the previous Section can be elaborated 
further in relation to Kolmogorov theory. Notice that the Kolmogorov interval in Fig. 
approximately coincides with the interval of wavenumbers in Fig. [^b) occupied by most of 
the pancake structures (except for the few largest local maximums with considerably larger 
wavenumbers). The Kolmogorov theory suggests the power-law (|(^u|) oc £^3 fQj- mean 
velocity variation at spatial scales i in the inertial interval, see e.g. [HI [I2|. Similarly, it 
implies (|5tc;|) oc for the vorticity variation. One can guess that the last relation is 

mainly determined by the small-scale high vorticity regions. In our simulations, that can 
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FIG. 12: (Color on-line) Relation between the vorticity local maximums Cc;max(^) and the respec¬ 
tive characteristic lengths li{t) during the evolution of the pancake structures: (a) for the global 
maximum, (b) for all local maximums. Red circles mark the local maximums at the final time 


—2/3 

t = 6.89. Dashed lines indicate the power-law Wmax oc ' in Eq. (17). 


be seen as describing the initial stage of tnrbulent dynamics, these high vorticity regions 
appear in the form of the pancake structures. As we argued in the previous Section, each 
pancake is characterized by the spatial scale ^i(t) oc (the other two scales ^2 and ^3 

remain close to 1 ) and the vorticity variation near the pancake is |hc<;| ~ oO m^^ U) oc 
Thus, we can test whether relation 

oc , (17) 


holds during the evolution of the pancake structures. 

(a), the global vorticity maximum indeed evolves along the 


As shown in Fig. 


12 


law. This can be checked additionally by the ratio of the characteristic times and 
that, according to relation 0 . should be equal to = 2/3. Our simulation provides 

the close relation Ti/T^^ 0.7 for the global vorticity maximum. It is also close for the 
simulation of the second initial condition I 2 , which yields Te/T^^ ~ 0.63, see Appendix [B| As 
one can see from Fig. [I^b), where the values of Umax are plotted versus the characteristic 
lengths ii for all the local maximums, most local maximums have the tendency to follow the 
power law 0 asymptotically. Note that the observed behavior agrees with the asymptotic 


pancake dynamics in accordance with Eqs. (13)-(14). However, this model allows any value 
for the ratio Ti/T^^. The blowup scenario based on the vortex lines breaking |26l[59j leads to 
the 2/3 ratio for the pancake structure Ti/T^^ and, thus, it can be considered as a possible 
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theoretical justification if extended to the case of the exponential vorticity growth. It is 
remarkable that local maximums follow approximately the same scaling law 0 , with some 
shared constant prefactor before ^ and also that all local maximums taken at fixed time 
are distributed around this law, see red points in Fig. [I^b). These two properties indicate 
strong correlation of pancakes in the process of energy transport to small scales at initial 
stages of turbulent flow. 


VI. CONCLUSIONS 


In this work we performed the systematic numerical study of high vorticity structures 
that develop in the 3D incompressible Euler equations from smooth initial conditions of 
finite energy. Being motivated by the open problem of the finite-time blowup and its role 
for the developed turbulence, we are led to two important observations. First, we show that 
the exponential growth of vorticity, which is typically observed within pancakes (thin and 
wide vortex sheets), is compatible with the formation of Kolmogorov turbulent spectra in 
the fully inviscid flow, i.e., before the viscous scales get excited. Second, we show that the 
pancake structures, which have self-similar dynamics and develop in increasing number, play 
the crucial role in formation of the energy cascade to small scales. 

We demonstrated that the thickness of pancake-like regions of high vorticity decreases 
exponentially in time, f'i(t) oc while the other two dimensions do not change consid¬ 
erably, £2 ~ ~ 1- At the same time the local vorticity maximum grows exponentially, 

<^max(^) oc During the evolution, the relation oc resembling the Kol¬ 
mogorov scaling law holds approximately in agreement with the ratio ~ 2/3 between 

the characteristic times of the pancake compression Tn and the vorticity growth Since 
the pancake evolution is governed by the two different time scales Tn and this behavior 
does not agree with the pancake model proposed in j35]. However, the modified model 


(13)-(14) satisfies the Euler equations for the leading terms and adequately describes the 


observed dynamics. The total number of pancake structures, estimated by the number of 
local vorticity maximums, increases with time. We demonstrate that at late times most of 
the pancakes are distributed densely across the corresponding interval of wavenumbers. 

We clearly observe the formation of the Kolmogorov energy spectrum oc 
the inviscid system, together with the exponential (i.e., no finite-time blowup) vorticity 
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growth. The interval with Kolmogorov scaling grows with time and extends to a decade 
of wavenumbers at the end of the simulations. The energy spectrum Ek{t) changes weakly 
with time in this region in contrast to vast changes at larger wavenumbers. 

Thin pancake structures in physical space generate strongly anisotropic vorticity field in 
Fourier space in the form of “jets”, which are extended in the directions perpendicular to 
the pancakes. Within these jets the Fourier components of the flow are large in comparison 
with the remaining background. We demonstrate that these jets occupy a small fraction 
of the entire spectral band, but provide the leading contribution to the energy spectrum 
of the system. This means that the energy transfer to small scales is performed through 
the evolution of the pancake structures, and that the Kolmogorov energy spectrum may be 
attributed to collective behavior of the pancakes. 
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Appendix A: Initial conditions 


TABLE I: Nonzero coefficients in Eq. (§ for the initial vorticity field Ii with the average energy 
density E/{27r)^ ^ 0.54 and the average helicity density i4/(27r)^ 1.05. Einal time for this 

simulation is t = 6.89 and the final grid is 486 x 1024 x 2048. The simulation is affected by aliasing 
starting from t 5.13. 


h 

Ah 

Bh 

(-1,0,2) 

(0.0065641, 0.0027931, 0.003282) 

(0.0044136, 0.0056271, 0.0022068) 

(0,0,0) 

(0.065101, 0.0005801, -0.064109) 

(0.0045744, -0.022895, 0.18392) 

(0,0,1) 

(0, 1, 0) 

(1, 0, 0) 

(0,0,2) 

(0, 0.01, 0) 

(0.01, 0, 0) 

(0,1,0) 

(0.21204, 0, -0.070625) 

(-0.14438, 0, 0.23298) 

(0,1,1) 

(0.045977, -0.010151, 0.010151) 

(0.041942, 0.040326, -0.040326) 

(0,2,0) 

(0.005, 0, 0) 

(0, 0, 0.005) 

(1,0,0) 

(0, 0, 0.1) 

(0, 0.1, 0) 

(1,0,1) 

(-0.046112, 0.017081, 0.046112) 

(-0.0097784, 0.020122, 0.0097784) 

(1,1,2) 

(-0.0034664, 0.0049556, -0.00074462) 

(-0.0059316, -0.0010472, 0.0034894) 

(2,0,0) 

(0, 0, 0.02) 

(0, 0.02, 0) 


TABLE II: Nonzero coefficients in Eq. Q for the initial vorticity field I 2 with the average energy 
density E/ (27r)^ 0.51 and the average helicity density fl/(27r)^ 1. Einal time for this simulation 

is t = 7.77 and the final grid is 1152 x 384 x 2304. The simulation is affected by aliasing starting 
from t ~ 6.25. 


h 

Ah 

Bh 

(-1,0,1) 

(-0.040618, 0.039651, -0.040618) 

(-0.030318, 0.064657, -0.030318) 

(0,0,0) 

(0.067751, -0.1311, -0.11256) 

(-0.082614, -0.0364, 0.18932) 

(0,0,1) 

(0, 1, 0) 

(1, 0, 0) 

(0,1,1) 

(0.0062549, 0.044315, -0.044315) 

(0.034983, -0.014521, 0.014521) 

(1,0,0) 

(0, 0.079395, 0.07027) 

(0, 0.099411, 0.012762) 

(1,1,0) 

(-0.047174, 0.047174, -0.045572) 

(-0.049622, 0.049622, 0.001773) 
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Appendix B: Simulation results for the initial condition I2 


In this Appendix we provide some results for the simulation with initial condition I 2 , 
which reaches the time t = 7.77 on the grid 1152 x 384 x 2304. The region of global vorticity 
maximum represents at final time a very thin pancake structure, as shown in Fig. [T^ by 
the numerically computed isosurface |tj(r)| = O.StUmax in the local coordinates ( 01 , 02 , 03 ), 
see Eq. a. The number of local maximums increases with time, from 2 at t = 0 to 13 
at the end of the simulation. The values of vorticity modulus at maximum points tend to 
grow exponentially in time, u)ma.xit) oc with different but relatively close values of the 

characteristic times Fig. [T^a). The thickness of the associated high vorticity regions 
decays nearly exponentially, £i{t) oc while the other two scales £2 and £3 remain the 

same or decrease slightly in time. Fig. [l4)(b). 

The gradual formation of the Kolmogorov region oc in the energy spectrum is 

shown in Fig. [I^a). This region extends to 5 ^ 15 at the end of the simulation, and 

corresponds to the “frozen” part of the spectrum where Ekit) changes slightly in time, in 
contrast to the vast changes at larger wavenumbers. The size of this interval is smaller than 
for the first simulation in Fig.|^ but the same tendency is clearly observed. Fig.[I^b) shows 
the characteristic wavenumbers ki = \/£\ for local maximums in decreasing order. The 
first three local maximums propagate much faster to higher wavenumbers, while the other 
local maximums fill densely the interval from large to medium scales, which approximately 
coincides with the Kolmogorov region in the energy spectrum. 

Fig. [TWa) shows the isosurface of renormalized Fourier components of vorticity |c5(k)| = 


0 . 2 , see Eq. (16), where the very thin interior part corresponds to larger vorticity. This 


isosurface is aligned close to the eigenvector wi computed at the global vorticity maximum. 
At smaller wavenumbers k ^ 15, corresponding to the Kolmogorov interval in Fig. [l^a), 
the isosurface represents a union of several jets, Fig. [T^b), that are clearly related to other 
local vorticity maximums (the figure shows the directions of the first, second and fourth 
largest local maximums, while the third local maximum yields the direction very close to 
that for the first one). The small interior part of the isosurface |t5(k)| = 0.2 dominates in 
the energy spectrum, Fig. 113 a). The contributions of the jets to the energy spectrum are 
comparable at sufficiently small wavenumbers from the Kolmogorov region, while at larger 
wavenumbers (both inside and outside the Kolmogorov region) the contribution from the 
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leading jet becomes dominant, Fig. [T7|(b). As shown in Fig. 18, vorticity maximums tend 
to evolve along the cjmax(^) power-law, with rather close values of the constant 

prefactors before for most of them. Thus, the simulation for the second initial condition 
leads to the same conclusions on the formation and structure of the Kolmogorov turbulent 
spectrum, as deduced for the first simulation in the main text of the paper. 



FIG. 13: (Color on-line) Isosurface of constant vorticity |a;| = O.ScJmax in the local coordinates 
(ai, a 2 , aa) at the final time of the simulation, t — 7.77. Note much smaller vertical scale. 



FIG. 14: (Color on-line) (a) Evolution of local vorticity maximums (logarithmic vertical scale). 
The dashed red line indicates the exponential slope oc with characteristic time = 2.7. 
(b) Evolution of characteristic spatial scales £i (black), £2 (blue) and £3 (red) of local maximums, 
see Eq. Q. The dashed red line indicates the exponential slope with characteristic time 


Te = 1.7. 
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FIG. 15: (Color on-line) (a) Energy spectrum E]^{t) at different times. Straight line above the 
curves indicates the slope of the Kolmogorov power-law, oc (b) Characteristic wavenum¬ 


bers ki estimated by Eq. (15) for the local maximums versus local maximum index numbers at 
different times. Local maximums are sorted in decreasing order of ki. Large circles mark the global 
maximums. 



FIG. 16: (Color on-line) Isosurface |d;(k)| = 0.2 of the normalized vorticity field in Fourier 
space at final time t — 7.77. (a) The jet is aligned with the eigenvector wi (dashed black line), 
which is the normal direction of the pancake structure at the global vorticity maximum in physical 
space, (b) The closer view with the eigenvector wi (solid black arrow) for the global vorticity 
maximum and the similar eigenvectors for the second (dashed blue arrow) and the fourth (dash- 
dot red arrow) largest local maximums. 
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FIG. 17: (Color on-line) Energy spectrum at the end of the simulation, (a) Calculated inside 
the isosurface |ci;(k)| > 0.2 (solid black), outside this isosurface (dash-dot blue) and the sum of the 
two (dashed red), (b) Calculated inside the isosurface |ci;(k)| > 0.2 (thick solid black), inside the 
first jet (dashed blue), inside the second and third jets (thin solid green), and the sum of the three 
jets (dash-dot red). Black line above the curves indicates the slope of the Kolmogorov power-law, 
Ek oc 



FIG. 18: (Color on-line) Relation between the vorticity local maximums Cc;max(^) and the respec¬ 
tive characteristic lengths li{t) during the evolution of the pancake structures: (a) for the global 
maximum, (b) for all local maximums. Red circles mark the local maximums at the final time 


— 2/3 

t = 7.77. Dashed lines indicate the power-law Cc;max oi ' in Eq. (17). 
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